R code (version 2026.01.2+418):

* Note: all paths used are absolute; in case of interest in repeat the analyses, the data path must be changed

### CONFIRMATORY FACTOR ANALYSIS CODE

# ========================================================================
# 1.CFA Analysis with Polychoric Correlations for Likert Data - 1 Factor Solution
# ========================================================================

# --------------------------
# Load Required Libraries
# --------------------------
library(haven)       # For reading SPSS files
library(lavaan)      # For CFA modeling
library(semPlot)     # For path diagrams
library(semTools)    # For additional SEM tools including reliability
library(psych)       # For polychoric correlations
library(dplyr)

# --------------------------
# Load and Prepare Data
# --------------------------
data_path <- "C:/Users/gigle/Downloads/GAS_data_2026.sav"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA1", "GASA2", "GASA3", "GASA4", "GASA5", "GASA6", "GASA7")]

# Convert all variables to numeric
gasa_data[] <- lapply(gasa_data, as.numeric)

# --------------------------
# Calculate Polychoric Correlations (appropriate for ordinal Likert data)
# --------------------------
poly_cor <- psych::polychoric(gasa_data)
cat("Polychoric Correlation Matrix:\n")
print(round(poly_cor$rho, 3))

# --------------------------
# Define and Run CFA Model using WLSMV
# --------------------------
cfa_model <- '
# Measurement model
PG =~ GASA1 + GASA2 + GASA3 + GASA4 + GASA5 + GASA6 + GASA7
'

cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares Mean and Variance adjusted (for ordinal data)
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# --------------------------
# Calculate All Required Fit Indices
# --------------------------
# Basic fit indices
fit_measures <- fitMeasures(cfa_result,
                            c("chisq", "df", "pvalue",
                              "cfi", "tli",
                              "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                              "srmr", "aic", "bic"))

# Additional indices
df_k <- fitMeasures(cfa_result, "df")            # DF of the fitted model
df_0 <- fitMeasures(cfa_result, "baseline.df")   # DF of the null model
cfi <- fitMeasures(cfa_result, "cfi")            # CFI
nfi <- fitMeasures(cfa_result, "nfi")            # Normed Fit Index

# GFI and AGFI if available with WLSMV
tryCatch({
  gfi <- fitMeasures(cfa_result, "gfi")
  agfi <- fitMeasures(cfa_result, "agfi")
}, error = function(e) {
  # If not available, calculate manually or use NA
  gfi <- NA
  agfi <- NA
  cat("Note: GFI and AGFI not directly available with WLSMV estimator\n")
})

# Calculate PCFI (Parsimony-adjusted CFI) - (df_model / df_null) * CFI
pcfi <- (df_k / df_0) * cfi

# Calculate Chi-square/df ratio
chisq_df_ratio <- fit_measures["chisq"] / fit_measures["df"]

# --------------------------
# Reliability and Validity Measures
# --------------------------
reliability_result <- semTools::reliability(cfa_result)

# --------------------------
# AVE Calculation
# --------------------------
ave_pg <- as.numeric(semTools::AVE(cfa_result)["PG"])

# --------------------------
# Calculate Modification Indices
# --------------------------
mod_indices <- modificationIndices(cfa_result, sort = TRUE, minimum.value = 3.84)

# --------------------------
# Parameter Estimates
# --------------------------
# R-squared values for indicators
r2_values <- lavInspect(cfa_result, "r2")

# --------------------------
# Create Summary Table
# --------------------------
# Create a data frame for the summary table
comparison_table <- data.frame(
  Index = c(
    "Chi-square", "df", "Chi-square/df",
    "CFI", "TLI", "NFI", "RMSEA", "RMSEA 90% CI", "SRMR",
    "GFI", "AGFI", "PCFI",
    "AIC", "BIC",
    "AVE PG"   # ← SIN COMA AQUÍ
  ),
  CFA = c(
    round(fit_measures["chisq"], 3),
    round(fit_measures["df"], 0),
    round(chisq_df_ratio, 3),
    round(fit_measures["cfi"], 3),
    round(fit_measures["tli"], 3),
    round(nfi, 3),
    round(fit_measures["rmsea"], 3),
    paste0("[", round(fit_measures["rmsea.ci.lower"], 3), ", ", round(fit_measures["rmsea.ci.upper"], 3), "]"),
    round(fit_measures["srmr"], 3),
    ifelse(is.na(gfi), NA, round(gfi, 3)),
    ifelse(is.na(agfi), NA, round(agfi, 3)),
    round(pcfi, 3),
    ifelse("aic" %in% names(fit_measures), round(fit_measures["aic"], 1), NA),
    ifelse("bic" %in% names(fit_measures), round(fit_measures["bic"], 1), NA),
    round(ave_pg, 3)
  )
)

# Print the final table
cat("\n--- Complete Fit Indices ---\n")
print(comparison_table, row.names = FALSE)

# Simple CFA Plot Code
# Load required libraries
library(lavaan)  # For CFA modeling
library(semPlot) # For path diagrams
library(haven)   # For reading SPSS data

# Data path
data_path <- "C:/Users/gigle/Documents/GASA_data"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA1", "GASA2", "GASA3", "GASA4", "GASA5", "GASA6", "GASA7")]

# Define the CFA model
cfa_model <- '
  # Measurement model
  PG =~ GASA1 + GASA2 + GASA3 + GASA4 + GASA5 + GASA6 + GASA7
'

# Run the CFA model - using WLSMV for ordinal data
cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares for ordinal data
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# Create the plot
semPaths(cfa_result,
         what = "std",           # Plot standardized solution
         layout = "tree2",       # Layout type
         edge.label.cex = 1.0,   # Size of edge labels
         whatLabels = "std",     # Show standardized values
         sizeMan = 8,            # Size of manifest variables
         sizeLat = 12,           # Size of latent variables
         curve = 1.5,            # Curve of the edges
         edge.color = "black",   # Color of edges
         border.width = 1.5,     # Width of borders
         fade = FALSE,           # Don't fade paths by strength
         title = FALSE,          # Don't include a title
         residuals = TRUE,       # Show residuals
         intercepts = FALSE)     # Don't show intercepts

# To save the plot as a PDF file
pdf("cfa_model_plot_onefactor.pdf", width = 10, height = 8)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "black",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = TRUE,
         intercepts = FALSE)
dev.off()

# To save the plot as a PNG file (higher resolution)
png("MODELO_CFA_GASA_1_FACTOR.png", width = 1200, height = 800, res = 120)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "blue",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = FALSE,
         intercepts = FALSE)
dev.off()

# Finally, print R² values
r2_values <- lavInspect(cfa_result, "r2")
print("R-squared values for each item:")
print(round(r2_values, 2))

# --------------------------
# Calculate Pseudo-AIC and BIC for WLSMV (approximate formulas  since AIC/BIC are not directly available for WLSMV)
# --------------------------
npar <- fitMeasures(cfa_result, "npar")          # Number of estimated parameters
nobs <- lavInspect(cfa_result, "ntotal")         # Sample size
pseudo_AIC <- fit_measures["chisq"] + 2 * npar   # Pseudo-AIC formula
pseudo_BIC <- fit_measures["chisq"] + npar * log(nobs)  # Pseudo-BIC formula
print(pseudo_AIC)
print(pseudo_BIC)


# ========================================================================
# 2.CFA Analysis with Polychoric Correlations for Likert Data - 2 Factor Solution
# ========================================================================

# --------------------------
# Load Required Libraries
# --------------------------
library(haven)       # For reading SPSS files
library(lavaan)      # For CFA modeling
library(semPlot)     # For path diagrams
library(semTools)    # For additional SEM tools including reliability
library(psych)       # For polychoric correlations
library(dplyr)

# --------------------------
# Load and Prepare Data
# --------------------------
data_path <- "C:/Users/gigle/Downloads/GAS_data_2026.sav"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA1", "GASA2", "GASA3", "GASA4", "GASA5",  "GASA6", "GASA7")]

# Convert all variables to numeric
gasa_data[] <- lapply(gasa_data, as.numeric)

# -------------------------------------------------------------------------
# Calculate Polychoric Correlations - (appropriate for ordinal Likert data)
# -------------------------------------------------------------------------
poly_cor <- psych::polychoric(gasa_data)
cat("Polychoric Correlation Matrix:\n")
print(round(poly_cor$rho, 3))

# --------------------------
# Define and Run CFA Model using WLSMV
# --------------------------
# Two-factor model
cfa_model <- '
# Measurement model
F1 =~ GASA1 + GASA2 + GASA3
F2 =~ GASA4 + GASA5 + GASA6 + GASA7
'

cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares Mean and Variance adjusted (for ordinal data)
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# --------------------------
# Calculate All Required Fit Indices
# --------------------------
# Basic fit indices
fit_measures <- fitMeasures(cfa_result,
                            c("chisq", "df", "pvalue",
                              "cfi", "tli",
                              "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                              "srmr", "aic", "bic"))

# Additional indices
df_k <- fitMeasures(cfa_result, "df")            # DF of the fitted model
df_0 <- fitMeasures(cfa_result, "baseline.df")   # DF of the null model
cfi <- fitMeasures(cfa_result, "cfi")            # CFI
nfi <- fitMeasures(cfa_result, "nfi")            # Normed Fit Index

# Try to get GFI and AGFI if available with WLSMV
tryCatch({
  gfi <- fitMeasures(cfa_result, "gfi")
  agfi <- fitMeasures(cfa_result, "agfi")
}, error = function(e) {
  # If not available, calculate manually or use NA
  gfi <- NA
  agfi <- NA
  cat("Note: GFI and AGFI not directly available with WLSMV estimator\n")
})

# Calculate PCFI (Parsimony-adjusted CFI)
pcfi <- (df_k / df_0) * cfi

# Calculate Chi-square/df ratio
chisq_df_ratio <- fit_measures["chisq"] / fit_measures["df"]

# --------------------------
# Reliability and Validity Measures
# --------------------------
# Calculate reliability measures
reliability_result <- semTools::reliability(cfa_result)

# --------------------------
# AVE Calculation for both factors
# --------------------------
ave_vals <- semTools::AVE(cfa_result)
ave_f1 <- as.numeric(ave_vals["F1"])
ave_f2 <- as.numeric(ave_vals["F2"])

# --------------------------
# Calculate HTMT (Heterotrait-Monotrait ratio)
# --------------------------
factor_corr <- lavInspect(cfa_result, "cor.lv")
f1_f2_corr <- factor_corr[1, 2]

htmt_mat <- semTools::htmt(
  cfa_model,
  data = gasa_data,
  ordered = names(gasa_data),
  missing = "listwise",
  htmt2 = TRUE
)
htmt_value <- as.numeric(htmt_mat[1, 2])

# --------------------------
# Calculate Modification Indices
# --------------------------
mod_indices <- modificationIndices(cfa_result, sort = TRUE, minimum.value = 3.84)

# --------------------------
# Parameter Estimates
# --------------------------
# R-squared values for indicators
r2_values <- lavInspect(cfa_result, "r2")

# --------------------------
# Summary of the model
# --------------------------
summary(cfa_result, fit.measures = TRUE, standardized = TRUE)

# --------------------------
# Create Summary Table
# --------------------------
# Create a data frame for the summary table
comparison_table <- data.frame(
  Index = c(
    "Chi-square", "df", "Chi-square/df",
    "CFI", "TLI", "NFI", "RMSEA", "RMSEA 90% CI", "SRMR",
    "GFI", "AGFI", "PCFI",
    "AIC", "BIC",
    "AVE F1", "AVE F2",
    "Factor Correlation (F1-F2)",
    "HTMT ratio"
  ),
  CFA = c(
    round(fit_measures["chisq"], 3),
    round(fit_measures["df"], 0),
    round(chisq_df_ratio, 3),
    round(fit_measures["cfi"], 3),
    round(fit_measures["tli"], 3),
    round(nfi, 3),
    round(fit_measures["rmsea"], 3),
    paste0("[", round(fit_measures["rmsea.ci.lower"], 3), ", ", round(fit_measures["rmsea.ci.upper"], 3), "]"),
    round(fit_measures["srmr"], 3),
    ifelse(is.na(gfi), NA, round(gfi, 3)),
    ifelse(is.na(agfi), NA, round(agfi, 3)),
    round(pcfi, 3),
    ifelse("aic" %in% names(fit_measures), round(fit_measures["aic"], 1), NA),
    ifelse("bic" %in% names(fit_measures), round(fit_measures["bic"], 1), NA),
    round(ave_f1, 3),
    round(ave_f2, 3),
    round(f1_f2_corr, 3),
    round(htmt_value, 3)
  )
)

cat("\n--- Complete Fit Indices ---\n")
print(comparison_table, row.names = FALSE)

# --------------------------
# Create CFA Path Diagram
# --------------------------
# Simple CFA Plot Code
# Load required libraries (already loaded above, but included here for completeness)
library(lavaan)  # For CFA modeling
library(semPlot) # For path diagrams
library(haven)   # For reading SPSS data

# Create the plot
semPaths(cfa_result,
         what = "std",           # Plot standardized solution
         layout = "tree2",       # Layout type
         edge.label.cex = 1.0,   # Size of edge labels
         whatLabels = "std",     # Show standardized values
         sizeMan = 8,            # Size of manifest variables
         sizeLat = 12,           # Size of latent variables
         curve = 1.5,            # Curve of the edges
         edge.color = "black",   # Color of edges
         border.width = 1.5,     # Width of borders
         fade = FALSE,           # Don't fade paths by strength
         title = FALSE,          # Don't include a title
         residuals = TRUE,       # Show residuals
         intercepts = FALSE)     # Don't show intercepts

# To save the plot as a PDF file
pdf("cfa_model_plot_twofactors.pdf", width = 10, height = 8)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "black",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = FALSE,
         intercepts = FALSE)
dev.off()

# To save the plot as a PNG file (higher resolution)
png("MODELO_CFA_GASA_2_FACTORS.png", width = 1200, height = 800, res = 120)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "purple",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = TRUE,
         intercepts = FALSE)
dev.off()

# Finally, print R² values
r2_values <- lavInspect(cfa_result, "r2")
print("R-squared values for each item:")
print(round(r2_values, 2))


# --------------------------
# Calculate Pseudo-AIC and BIC for WLSMV
# --------------------------
npar <- fitMeasures(cfa_result, "npar")          # Number of estimated parameters
nobs <- lavInspect(cfa_result, "ntotal")         # Sample size
pseudo_AIC <- fit_measures["chisq"] + 2 * npar   # Pseudo-AIC formula
pseudo_BIC <- fit_measures["chisq"] + npar * log(nobs)  # Pseudo-BIC formula
print(pseudo_AIC)
print(pseudo_BIC)


# ========================================================================
# 3.CFA Analysis with Polychoric Correlations for Likert Data - CORE model
# ========================================================================

# --------------------------
# Load Required Libraries
# --------------------------
library(haven)       # For reading SPSS files
library(lavaan)      # For CFA modeling
library(semPlot)     # For path diagrams
library(semTools)    # For additional SEM tools including reliability
library(psych)       # For polychoric correlations
library(dplyr)

# --------------------------
# Load and Prepare Data
# --------------------------
data_path <- "C:/Users/gigle/Downloads/GAS_data_2026.sav"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA4", "GASA5",  "GASA6", "GASA7")]

# Convert all variables to numeric
gasa_data[] <- lapply(gasa_data, as.numeric)

# --------------------------
# Calculate Polychoric Correlations
# --------------------------
# Compute polychoric correlation matrix (appropriate for ordinal Likert data)
poly_cor <- psych::polychoric(gasa_data)
cat("Polychoric Correlation Matrix:\n")
print(round(poly_cor$rho, 3))

# --------------------------
# Define and Run CFA Model using WLSMV
# --------------------------
cfa_model <- '
# Measurement model
PG =~ GASA4 + GASA5 + GASA6 + GASA7
'

cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares Mean and Variance adjusted (for ordinal data)
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# --------------------------
# Calculate All Required Fit Indices
# --------------------------
# Basic fit indices
fit_measures <- fitMeasures(cfa_result,
                            c("chisq", "df", "pvalue",
                              "cfi", "tli",
                              "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                              "srmr", "aic", "bic"))

# Additional indices
df_k <- fitMeasures(cfa_result, "df")            # DF of the fitted model
df_0 <- fitMeasures(cfa_result, "baseline.df")   # DF of the null model
cfi <- fitMeasures(cfa_result, "cfi")            # CFI
nfi <- fitMeasures(cfa_result, "nfi")            # Normed Fit Index

# Try to get GFI and AGFI if available with WLSMV
tryCatch({
  gfi <- fitMeasures(cfa_result, "gfi")
  agfi <- fitMeasures(cfa_result, "agfi")
}, error = function(e) {
  # If not available, calculate manually or use NA
  gfi <- NA
  agfi <- NA
  cat("Note: GFI and AGFI not directly available with WLSMV estimator\n")
})

# Calculate PCFI (Parsimony-adjusted CFI)
# Correct formula: PCFI = (df_model / df_null) * CFI
pcfi <- (df_k / df_0) * cfi

# Calculate Chi-square/df ratio
chisq_df_ratio <- fit_measures["chisq"] / fit_measures["df"]

# --------------------------
# Reliability and Validity Measures
# --------------------------
# Calculate reliability measures
reliability_result <- semTools::reliability(cfa_result)

# --------------------------
# AVE Calculation
# --------------------------
ave_pg <- as.numeric(semTools::AVE(cfa_result)["PG"])

# --------------------------
# Calculate Modification Indices
# --------------------------
mod_indices <- modificationIndices(cfa_result, sort = TRUE, minimum.value = 3.84)

# --------------------------
# Parameter Estimates
# --------------------------
# R-squared values for indicators
r2_values <- lavInspect(cfa_result, "r2")

# --------------------------
# Create Summary Table
# --------------------------
# Create a data frame for the summary table
comparison_table <- data.frame(
  Index = c(
    "Chi-square", "df", "Chi-square/df",
    "CFI", "TLI", "NFI", "RMSEA", "RMSEA 90% CI", "SRMR",
    "GFI", "AGFI", "PCFI",
    "AIC", "BIC",
    "AVE PG"
  ),
  CFA = c(
    round(fit_measures["chisq"], 3),
    round(fit_measures["df"], 0),
    round(chisq_df_ratio, 3),
    round(fit_measures["cfi"], 3),
    round(fit_measures["tli"], 3),
    round(nfi, 3),
    round(fit_measures["rmsea"], 3),
    paste0("[", round(fit_measures["rmsea.ci.lower"], 3), ", ", round(fit_measures["rmsea.ci.upper"], 3), "]"),
    round(fit_measures["srmr"], 3),
    ifelse(is.na(gfi), NA, round(gfi, 3)),
    ifelse(is.na(agfi), NA, round(agfi, 3)),
    round(pcfi, 3),
    ifelse("aic" %in% names(fit_measures), round(fit_measures["aic"], 1), NA),
    ifelse("bic" %in% names(fit_measures), round(fit_measures["bic"], 1), NA),
    round(ave_pg, 3)
  )
)

cat("\n--- Complete Fit Indices ---\n")
print(comparison_table, row.names = FALSE)

# Simple CFA Plot Code
# Load required libraries
library(lavaan)  # For CFA modeling
library(semPlot) # For path diagrams
library(haven)   # For reading SPSS data

# Specify your data path - again, just in case
data_path <- "C:/Users/gigle/Documents/GASA_data"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA4", "GASA5",  "GASA6", "GASA7")]

# Define the CFA model
cfa_model <- '
  # Measurement model
  PG =~ GASA4 + GASA5 + GASA6 + GASA7
'

# Run the CFA model - using WLSMV for ordinal data
cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares for ordinal data
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# Create the plot
semPaths(cfa_result,
         what = "std",           # Plot standardized solution
         layout = "tree2",       # Layout type
         edge.label.cex = 1.0,   # Size of edge labels
         whatLabels = "std",     # Show standardized values
         sizeMan = 8,            # Size of manifest variables
         sizeLat = 12,           # Size of latent variables
         curve = 1.5,            # Curve of the edges
         edge.color = "black",   # Color of edges
         border.width = 1.5,     # Width of borders
         fade = FALSE,           # Don't fade paths by strength
         title = FALSE,          # Don't include a title
         residuals = TRUE,       # Show residuals
         intercepts = FALSE)     # Don't show intercepts

# To save the plot as a PDF file
pdf("cfa_model_plot_onefactor.pdf", width = 10, height = 8)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "black",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = FALSE,
         intercepts = FALSE)
dev.off()

# To save the plot as a PNG file (higher resolution)
png("MODELO_CFA_GASA_1_FACTOR.png", width = 1200, height = 800, res = 120)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "blue",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = TRUE,
         intercepts = FALSE)
dev.off()

# Finally, print R² values
r2_values <- lavInspect(cfa_result, "r2")
print("R-squared values for each item:")
print(round(r2_values, 2))

### PSEUDO AIC and BIC

# --------------------------
# Calculate Pseudo-AIC and BIC for WLSMV
# --------------------------
npar <- fitMeasures(cfa_result, "npar")          # Number of estimated parameters
nobs <- lavInspect(cfa_result, "ntotal")         # Sample size
pseudo_AIC <- fit_measures["chisq"] + 2 * npar   # Pseudo-AIC formula
pseudo_BIC <- fit_measures["chisq"] + npar * log(nobs)  # Pseudo-BIC formula
print(pseudo_AIC)
print(pseudo_BIC)


# ========================================================================
# 4a.CFA Analysis with Polychoric Correlations for Likert Data - ICD-11 model
### -> As a model with 3 items and 1 factor has 0 degrees of freedom, so goodness of fit
# indices cannot be computed as it is just-identified: those indices would be artificially # perfect. The following code correspons to a 1-factor model used to calculate factor
# loadings and R2 values
# ========================================================================

# --------------------------
# Load Required Libraries
# --------------------------
library(haven)       # For reading SPSS files
library(lavaan)      # For CFA modeling
library(semPlot)     # For path diagrams
library(semTools)    # For additional SEM tools including reliability
library(psych)       # For polychoric correlations
library(dplyr)

# --------------------------
# Load and Prepare Data
# --------------------------
data_path <- "C:/Users/gigle/Downloads/GAS_data_2026.sav"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA4", "GASA6", "GASA7")]

# Convert all variables to numeric
gasa_data[] <- lapply(gasa_data, as.numeric)

# --------------------------
# Calculate Polychoric Correlations (appropriate for ordinal Likert data)
# --------------------------
poly_cor <- psych::polychoric(gasa_data)
cat("Polychoric Correlation Matrix:\n")
print(round(poly_cor$rho, 3))

# --------------------------
# Define and Run CFA Model using WLSMV
# --------------------------
cfa_model <- '
# Measurement model
PG =~ GASA4 + GASA6 + GASA7
'

cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares Mean and Variance adjusted (for ordinal data)
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# --------------------------
# Calculate All Required Fit Indices
# --------------------------
# Basic fit indices
fit_measures <- fitMeasures(cfa_result,
                            c("chisq", "df", "pvalue",
                              "cfi", "tli",
                              "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                              "srmr", "aic", "bic"))

# Additional indices
df_k <- fitMeasures(cfa_result, "df")            # DF of the fitted model
df_0 <- fitMeasures(cfa_result, "baseline.df")   # DF of the null model
cfi <- fitMeasures(cfa_result, "cfi")            # CFI
nfi <- fitMeasures(cfa_result, "nfi")            # Normed Fit Index

# GFI and AGFI if available with WLSMV
tryCatch({
  gfi <- fitMeasures(cfa_result, "gfi")
  agfi <- fitMeasures(cfa_result, "agfi")
}, error = function(e) {
  # If not available, calculate manually or use NA
  gfi <- NA
  agfi <- NA
  cat("Note: GFI and AGFI not directly available with WLSMV estimator\n")
})

# Calculate PCFI (Parsimony-adjusted CFI) - (df_model / df_null) * CFI
pcfi <- (df_k / df_0) * cfi

# Calculate Chi-square/df ratio
chisq_df_ratio <- fit_measures["chisq"] / fit_measures["df"]

# --------------------------
# Reliability and Validity Measures
# --------------------------
reliability_result <- semTools::reliability(cfa_result)

# --------------------------
# AVE Calculation
# --------------------------
ave_pg <- as.numeric(semTools::AVE(cfa_result)["PG"])

# --------------------------
# Calculate Modification Indices
# --------------------------
mod_indices <- modificationIndices(cfa_result, sort = TRUE, minimum.value = 3.84)

# --------------------------
# Parameter Estimates
# --------------------------
# R-squared values for indicators
r2_values <- lavInspect(cfa_result, "r2")

# --------------------------
# Create Summary Table
# --------------------------
# Create a data frame for the summary table
comparison_table <- data.frame(
  Index = c(
    "Chi-square", "df", "Chi-square/df",
    "CFI", "TLI", "NFI", "RMSEA", "RMSEA 90% CI", "SRMR",
    "GFI", "AGFI", "PCFI",
    "AIC", "BIC",
    "AVE PG"   # ← SIN COMA AQUÍ
  ),
  CFA = c(
    round(fit_measures["chisq"], 3),
    round(fit_measures["df"], 0),
    round(chisq_df_ratio, 3),
    round(fit_measures["cfi"], 3),
    round(fit_measures["tli"], 3),
    round(nfi, 3),
    round(fit_measures["rmsea"], 3),
    paste0("[", round(fit_measures["rmsea.ci.lower"], 3), ", ", round(fit_measures["rmsea.ci.upper"], 3), "]"),
    round(fit_measures["srmr"], 3),
    ifelse(is.na(gfi), NA, round(gfi, 3)),
    ifelse(is.na(agfi), NA, round(agfi, 3)),
    round(pcfi, 3),
    ifelse("aic" %in% names(fit_measures), round(fit_measures["aic"], 1), NA),
    ifelse("bic" %in% names(fit_measures), round(fit_measures["bic"], 1), NA),
    round(ave_pg, 3)
  )
)

# Print the final table
cat("\n--- Complete Fit Indices ---\n")
print(comparison_table, row.names = FALSE)

# Simple CFA Plot Code
# Load required libraries
library(lavaan)  # For CFA modeling
library(semPlot) # For path diagrams
library(haven)   # For reading SPSS data

# Data path - replace with your actual file path
data_path <- "C:/Users/gigle/Documents/GASA_data"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA4", "GASA6", "GASA7")]

# Define the CFA model
cfa_model <- '
  # Measurement model
  PG =~ GASA4 + GASA6 + GASA7
'

# Run the CFA model - using WLSMV for ordinal data
cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares for ordinal data
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# Create the plot
semPaths(cfa_result,
         what = "std",           # Plot standardized solution
         layout = "tree2",       # Layout type
         edge.label.cex = 1.0,   # Size of edge labels
         whatLabels = "std",     # Show standardized values
         sizeMan = 8,            # Size of manifest variables
         sizeLat = 12,           # Size of latent variables
         curve = 1.5,            # Curve of the edges
         edge.color = "black",   # Color of edges
         border.width = 1.5,     # Width of borders
         fade = FALSE,           # Don't fade paths by strength
         title = FALSE,          # Don't include a title
         residuals = TRUE,       # Show residuals
         intercepts = FALSE)     # Don't show intercepts

# To save the plot as a PDF file
pdf("cfa_model_plot_onefactor.pdf", width = 10, height = 8)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "black",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = TRUE,
         intercepts = FALSE)
dev.off()

# To save the plot as a PNG file (higher resolution)
png("MODELO_CFA_GASA_1_FACTOR.png", width = 1200, height = 800, res = 120)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "blue",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = FALSE,
         intercepts = FALSE)
dev.off()

# Finally, print R² values
r2_values <- lavInspect(cfa_result, "r2")
print("R-squared values for each item:")
print(round(r2_values, 2))

# --------------------------
# Calculate Pseudo-AIC and BIC for WLSMV (approximate formulas  since AIC/BIC are not directly available for WLSMV)
# --------------------------
npar <- fitMeasures(cfa_result, "npar")          # Number of estimated parameters
nobs <- lavInspect(cfa_result, "ntotal")         # Sample size
pseudo_AIC <- fit_measures["chisq"] + 2 * npar   # Pseudo-AIC formula
pseudo_BIC <- fit_measures["chisq"] + npar * log(nobs)  # Pseudo-BIC formula
print(pseudo_AIC)
print(pseudo_BIC)

# ========================================================================
# 4b.CFA Analysis with Polychoric Correlations for Likert Data - ICD-11 model:
### -> As a model with 3 items and 1 factor has 0 degrees of freedom, so goodness of fit
# indices cannot be computed as it is just-identified: those indices would be artificially # perfect. So, these results are approximate, calculated using a two-factor model instead, # with items 4,6 and 7 loading in one factor, and 1, 2, 3 and 5 in the other. The following
# code correspons to a 2-factor model used to calculate goodness of fit values, in order to maximize approximate comparability with the other factor solutions
# ========================================================================

# --------------------------
# Load Required Libraries
# --------------------------
library(haven)       # For reading SPSS files
library(lavaan)      # For CFA modeling
library(semPlot)     # For path diagrams
library(semTools)    # For additional SEM tools including reliability
library(psych)       # For polychoric correlations
library(dplyr)

# --------------------------
# Load and Prepare Data
# --------------------------
data_path <- "C:/Users/gigle/Documents/GASA_data"
data <- read_sav(data_path)

# Select only the GASA items
gasa_data <- data[, c("GASA1", "GASA2", "GASA3", "GASA4", "GASA5",  "GASA6", "GASA7")]

# Convert all variables to numeric
gasa_data[] <- lapply(gasa_data, as.numeric)

# -------------------------------------------------------------------------
# Calculate Polychoric Correlations - (appropriate for ordinal Likert data)
# -------------------------------------------------------------------------
poly_cor <- psych::polychoric(gasa_data)
cat("Polychoric Correlation Matrix:\n")
print(round(poly_cor$rho, 3))

# --------------------------
# Define and Run CFA Model using WLSMV
# --------------------------
# Two-factor model
cfa_model <- '
# Measurement model
F1 =~ GASA1 + GASA2 + GASA3 + GASA5
F2 =~ GASA4 + GASA6 + GASA7
'

cfa_result <- cfa(
  model = cfa_model,
  data = gasa_data,
  estimator = "WLSMV",      # Weighted Least Squares Mean and Variance adjusted (for ordinal data)
  ordered = names(gasa_data) # Specify all variables as ordinal
)

# --------------------------
# Calculate All Required Fit Indices
# --------------------------
# Basic fit indices
fit_measures <- fitMeasures(cfa_result,
                            c("chisq", "df", "pvalue",
                              "cfi", "tli",
                              "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                              "srmr", "aic", "bic"))

# Additional indices
df_k <- fitMeasures(cfa_result, "df")            # DF of the fitted model
df_0 <- fitMeasures(cfa_result, "baseline.df")   # DF of the null model
cfi <- fitMeasures(cfa_result, "cfi")            # CFI
nfi <- fitMeasures(cfa_result, "nfi")            # Normed Fit Index

# Try to get GFI and AGFI if available with WLSMV
tryCatch({
  gfi <- fitMeasures(cfa_result, "gfi")
  agfi <- fitMeasures(cfa_result, "agfi")
}, error = function(e) {
  # If not available, calculate manually or use NA
  gfi <- NA
  agfi <- NA
  cat("Note: GFI and AGFI not directly available with WLSMV estimator\n")
})

# Calculate PCFI (Parsimony-adjusted CFI)
pcfi <- (df_k / df_0) * cfi

# Calculate Chi-square/df ratio
chisq_df_ratio <- fit_measures["chisq"] / fit_measures["df"]

# --------------------------
# Reliability and Validity Measures
# --------------------------
# Calculate reliability measures
reliability_result <- semTools::reliability(cfa_result)

# --------------------------
# AVE Calculation for both factors
# --------------------------
ave_values <- semTools::AVE(cfa_result)

ave_f1 <- as.numeric(ave_values["F1"])
ave_f2 <- as.numeric(ave_values["F2"])

# --------------------------
# Calculate HTMT (Heterotrait-Monotrait ratio)
# --------------------------
# Get factor correlation
factor_corr <- lavInspect(cfa_result, "cor.lv")
f1_f2_corr <- factor_corr[1, 2]

# Calculate average monotrait correlations
f1_items <- c("GASA1", "GASA2", "GASA3", "GASA5")
f2_items <- c("GASA4", "GASA6", "GASA7")

# Extract the relevant correlations from the polychoric correlation matrix
f1_correlations <- poly_cor$rho[f1_items, f1_items]
f2_correlations <- poly_cor$rho[f2_items, f2_items]
hetero_correlations <- poly_cor$rho[f1_items, f2_items]

# Calculate average heterotrait correlations (excluding diagonals for monotrait)
avg_monotrait_f1 <- mean(f1_correlations[lower.tri(f1_correlations, diag = FALSE)])
avg_monotrait_f2 <- mean(f2_correlations[lower.tri(f2_correlations, diag = FALSE)])
avg_heterotrait <- mean(hetero_correlations)

# Calculate HTMT
htmt_value <- avg_heterotrait / sqrt(avg_monotrait_f1 * avg_monotrait_f2)

# --------------------------
# Calculate Modification Indices
# --------------------------
mod_indices <- modificationIndices(cfa_result, sort = TRUE, minimum.value = 3.84)

# --------------------------
# Parameter Estimates
# --------------------------
# R-squared values for indicators
r2_values <- lavInspect(cfa_result, "r2")

# --------------------------
# Summary of the model
# --------------------------
summary(cfa_result, fit.measures = TRUE, standardized = TRUE)

# --------------------------
# Create Summary Table
# --------------------------
# --------------------------
# Create Summary Table
# --------------------------
comparison_table <- data.frame(
  Index = c(
    "Chi-square", "df", "Chi-square/df",
    "CFI", "TLI", "NFI", "RMSEA", "RMSEA 90% CI", "SRMR",
    "GFI", "AGFI", "PCFI",
    "AIC", "BIC",
    "AVE F1", "AVE F2",
    "Factor Correlation (F1-F2)",
    "HTMT ratio"
  ),
  CFA = c(
    round(fit_measures["chisq"], 3),
    round(fit_measures["df"], 0),
    round(chisq_df_ratio, 3),
    round(fit_measures["cfi"], 3),
    round(fit_measures["tli"], 3),
    round(nfi, 3),
    round(fit_measures["rmsea"], 3),
    paste0("[", round(fit_measures["rmsea.ci.lower"], 3), ", ", round(fit_measures["rmsea.ci.upper"], 3), "]"),
    round(fit_measures["srmr"], 3),
    ifelse(is.na(gfi), NA, round(gfi, 3)),
    ifelse(is.na(agfi), NA, round(agfi, 3)),
    round(pcfi, 3),
    ifelse("aic" %in% names(fit_measures), round(fit_measures["aic"], 1), NA),
    ifelse("bic" %in% names(fit_measures), round(fit_measures["bic"], 1), NA),
    round(ave_f1, 3),
    round(ave_f2, 3),
    round(f1_f2_corr, 3),
    round(htmt_value, 3)
  )
)

cat("\n--- Complete Fit Indices ---\n")
print(comparison_table, row.names = FALSE)

# --------------------------
# Create CFA Path Diagram
# --------------------------
# Simple CFA Plot Code
# Load required libraries (already loaded above, but included here for completeness)
library(lavaan)  # For CFA modeling
library(semPlot) # For path diagrams
library(haven)   # For reading SPSS data

# Create the plot
semPaths(cfa_result,
         what = "std",           # Plot standardized solution
         layout = "tree2",       # Layout type
         edge.label.cex = 1.0,   # Size of edge labels
         whatLabels = "std",     # Show standardized values
         sizeMan = 8,            # Size of manifest variables
         sizeLat = 12,           # Size of latent variables
         curve = 1.5,            # Curve of the edges
         edge.color = "black",   # Color of edges
         border.width = 1.5,     # Width of borders
         fade = FALSE,           # Don't fade paths by strength
         title = FALSE,          # Don't include a title
         residuals = TRUE,       # Show residuals
         intercepts = FALSE)     # Don't show intercepts

# To save the plot as a PDF file
pdf("cfa_model_plot_twofactors.pdf", width = 10, height = 8)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "blue",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = FALSE,
         intercepts = FALSE)
dev.off()

# To save the plot as a PNG file (higher resolution)
png("MODELO_CFA_GASA_2_FACTORS.png", width = 1200, height = 800, res = 120)
semPaths(cfa_result,
         what = "std",
         layout = "tree2",
         edge.label.cex = 1.0,
         whatLabels = "std",
         sizeMan = 8,
         sizeLat = 12,
         curve = 1.5,
         edge.color = "purple",
         border.width = 1.5,
         fade = FALSE,
         title = FALSE,
         residuals = TRUE,
         intercepts = FALSE)
dev.off()

# Finally, print R² values
r2_values <- lavInspect(cfa_result, "r2")
print("R-squared values for each item:")
print(round(r2_values, 2))


# --------------------------
# Calculate Pseudo-AIC and BIC for WLSMV
# --------------------------
npar <- fitMeasures(cfa_result, "npar")          # Number of estimated parameters
nobs <- lavInspect(cfa_result, "ntotal")         # Sample size
pseudo_AIC <- fit_measures["chisq"] + 2 * npar   # Pseudo-AIC formula
pseudo_BIC <- fit_measures["chisq"] + npar * log(nobs)  # Pseudo-BIC formula
print(pseudo_AIC)
print(pseudo_BIC)


### OMEGA CALCULATION CODE

# =============================================================================
# RELIABILITY AND AVERAGE VARIANCE EXTRACTED 1 Factor Solution
# =============================================================================

library(lavaan)
library(boot)
library(haven)

# --------------------------
# Load data
# --------------------------
data_path <- "C:/Users/gigle/Documents/GASA_data"
datos <- read_sav(data_path)

items_all <- c("GASA1", "GASA2", "GASA3", "GASA4", "GASA5", "GASA6", "GASA7")

# --------------------------
# CFA model
# --------------------------
cfa_model <- '
  F1 =~ GASA1 + GASA2 + GASA3 + GASA4 + GASA5 + GASA6 + GASA7
'

# --------------------------
# FUNCTION: Omega from CFA
# --------------------------
omega_from_fit <- function(fit) {
 
  # standardized solution
  std <- standardizedSolution(fit)
 
  # loadings
  lambda <- std$est.std[std$op == "=~"]
 
  # residual variances
  theta <- std$est.std[std$op == "~~" & std$lhs %in% items_all & std$lhs == std$rhs]
 
  # omega formula
  omega <- (sum(lambda))^2 / ((sum(lambda))^2 + sum(theta))
 
  return(omega)
}

# --------------------------
# Bootstrap function
# --------------------------
boot_fun <- function(data, indices) {
 
  d <- data[indices, ]
 
  fit <- tryCatch(
    cfa(cfa_model,
        data = d,
        ordered = items_all,
        estimator = "WLSMV",
        std.lv = TRUE),
    error = function(e) NULL
  )
 
  if (is.null(fit)) return(NA)
 
  if (!lavInspect(fit, "converged")) return(NA)
 
  return(omega_from_fit(fit))
}

# --------------------------
# POINT ESTIMATE
# --------------------------
fit_point <- cfa(cfa_model,
                 data = datos,
                 ordered = items_all,
                 estimator = "WLSMV",
                 std.lv = TRUE)

omega_point <- omega_from_fit(fit_point)

cat("\nOmega:", round(omega_point, 3))

# --------------------------
# BOOTSTRAP
# --------------------------
set.seed(123)

boot_res <- boot(
  data = datos,
  statistic = boot_fun,
  R = 500,
  parallel = "multicore"
)

# clean NAs
boot_vals <- boot_res$t[!is.na(boot_res$t)]

# CI
ci <- quantile(boot_vals, c(0.025, 0.975))

cat("\n\nOmega with 95% CI:\n")
cat(sprintf("%.3f [%.3f, %.3f]", omega_point, ci[1], ci[2]))



# =============================================================================
# RELIABILITY AND AVERAGE VARIANCE EXTRACTED 2 Factor Solution
# =============================================================================

library(lavaan)
library(boot)
library(haven)

# --------------------------
# Load data
# --------------------------
data_path <- "C:/Users/gigle/Downloads/GAS_data_2026.sav"
datos <- read_sav(data_path)

items_F1 <- c("GASA1", "GASA2", "GASA3")
items_F2 <- c("GASA4", "GASA5", "GASA6", "GASA7")

items_all <- c(items_F1, items_F2)

# --------------------------
# CFA model
# --------------------------
cfa_model <- '
  F1 =~ GASA1 + GASA2 + GASA3
  F2 =~ GASA4 + GASA5 + GASA6 + GASA7
'

# --------------------------
# FUNCTION: Omega per factor
# --------------------------
omega_from_fit <- function(fit) {
 
  std <- standardizedSolution(fit)
 
  # -------- F1 --------
  lambda_F1 <- std$est.std[std$op == "=~" & std$lhs == "F1"]
  theta_F1  <- std$est.std[std$op == "~~" & std$lhs %in% items_F1 & std$lhs == std$rhs]
 
  omega_F1 <- (sum(lambda_F1))^2 / ((sum(lambda_F1))^2 + sum(theta_F1))
 
  # -------- F2 --------
  lambda_F2 <- std$est.std[std$op == "=~" & std$lhs == "F2"]
  theta_F2  <- std$est.std[std$op == "~~" & std$lhs %in% items_F2 & std$lhs == std$rhs]
 
  omega_F2 <- (sum(lambda_F2))^2 / ((sum(lambda_F2))^2 + sum(theta_F2))
 
  return(c(F1 = omega_F1, F2 = omega_F2))
}

# --------------------------
# Bootstrap function
# --------------------------
boot_fun <- function(data, indices) {
 
  d <- data[indices, ]
 
  fit <- tryCatch(
    cfa(cfa_model,
        data = d,
        ordered = items_all,
        estimator = "WLSMV",
        std.lv = TRUE),
    error = function(e) NULL
  )
 
  if (is.null(fit)) return(c(NA, NA))
  if (!lavInspect(fit, "converged")) return(c(NA, NA))
 
  return(omega_from_fit(fit))
}

# --------------------------
# POINT ESTIMATES
# --------------------------
fit_point <- cfa(cfa_model,
                 data = datos,
                 ordered = items_all,
                 estimator = "WLSMV",
                 std.lv = TRUE)

omega_point <- omega_from_fit(fit_point)

cat("\n===== Omega Point Estimates =====\n")
cat("\nF1:", round(omega_point["F1"], 3))
cat("\nF2:", round(omega_point["F2"], 3))

# --------------------------
# BOOTSTRAP
# --------------------------
set.seed(123)

boot_res <- boot(
  data = datos,
  statistic = boot_fun,
  R = 500,
  parallel = "multicore"
)

# limpiar NAs
boot_F1 <- boot_res$t[,1]
boot_F2 <- boot_res$t[,2]

boot_F1 <- boot_F1[!is.na(boot_F1)]
boot_F2 <- boot_F2[!is.na(boot_F2)]

# CIs
ci_F1 <- quantile(boot_F1, c(0.025, 0.975))
ci_F2 <- quantile(boot_F2, c(0.025, 0.975))

# --------------------------
# FINAL OUTPUT
# --------------------------
cat("\n\n===== Omega with 95% CI =====\n")

cat("\nF1:",
    sprintf("%.3f [%.3f, %.3f]",
            omega_point["F1"], ci_F1[1], ci_F1[2]))

cat("\nF2:",
    sprintf("%.3f [%.3f, %.3f]",
            omega_point["F2"], ci_F2[1], ci_F2[2]))


# =============================================================================
# RELIABILITY AND AVERAGE VARIANCE EXTRACTED CORE model
# =============================================================================

library(lavaan)
library(boot)
library(haven)

# --------------------------
# Load data
# --------------------------
data_path <- "C:/Users/gigle/Documents/GASA_data"
datos <- read_sav(data_path)

items_all <- c("GASA4", "GASA5", "GASA6", "GASA7")

# --------------------------
# CFA model
# --------------------------
cfa_model <- '
  F1 =~ GASA4 + GASA5 + GASA6 + GASA7
'

# --------------------------
# FUNCTION: Omega from CFA
# --------------------------
omega_from_fit <- function(fit) {
 
  # standardized solution
  std <- standardizedSolution(fit)
 
  # loadings
  lambda <- std$est.std[std$op == "=~"]
 
  # residual variances
  theta <- std$est.std[std$op == "~~" & std$lhs %in% items_all & std$lhs == std$rhs]
 
  # omega formula
  omega <- (sum(lambda))^2 / ((sum(lambda))^2 + sum(theta))
 
  return(omega)
}

# --------------------------
# Bootstrap function
# --------------------------
boot_fun <- function(data, indices) {
 
  d <- data[indices, ]
 
  fit <- tryCatch(
    cfa(cfa_model,
        data = d,
        ordered = items_all,
        estimator = "WLSMV",
        std.lv = TRUE),
    error = function(e) NULL
  )
 
  if (is.null(fit)) return(NA)
 
  if (!lavInspect(fit, "converged")) return(NA)
 
  return(omega_from_fit(fit))
}

# --------------------------
# POINT ESTIMATE
# --------------------------
fit_point <- cfa(cfa_model,
                 data = datos,
                 ordered = items_all,
                 estimator = "WLSMV",
                 std.lv = TRUE)

omega_point <- omega_from_fit(fit_point)

cat("\nOmega:", round(omega_point, 3))

# --------------------------
# BOOTSTRAP
# --------------------------
set.seed(123)

boot_res <- boot(
  data = datos,
  statistic = boot_fun,
  R = 500,
  parallel = "multicore"
)

# clean NAs
boot_vals <- boot_res$t[!is.na(boot_res$t)]

# CI
ci <- quantile(boot_vals, c(0.025, 0.975))

cat("\n\nOmega with 95% CI:\n")
cat(sprintf("%.3f [%.3f, %.3f]", omega_point, ci[1], ci[2]))


# =============================================================================
# RELIABILITY AND AVERAGE VARIANCE EXTRACTED ICD-11 model
# =============================================================================

library(lavaan)
library(boot)
library(haven)

# --------------------------
# Load data
# --------------------------
data_path <- "C:/Users/gigle/Documents/GASA_data"
datos <- read_sav(data_path)

items_all <- c("GASA4", "GASA6", "GASA7")

# --------------------------
# CFA model
# --------------------------
cfa_model <- '
  F1 =~ GASA4 + GASA6 + GASA7
'

# --------------------------
# FUNCTION: Omega from CFA
# --------------------------
omega_from_fit <- function(fit) {
 
  # standardized solution
  std <- standardizedSolution(fit)
 
  # loadings
  lambda <- std$est.std[std$op == "=~"]
 
  # residual variances
  theta <- std$est.std[std$op == "~~" & std$lhs %in% items_all & std$lhs == std$rhs]
 
  # omega formula
  omega <- (sum(lambda))^2 / ((sum(lambda))^2 + sum(theta))
 
  return(omega)
}

# --------------------------
# Bootstrap function
# --------------------------
boot_fun <- function(data, indices) {
 
  d <- data[indices, ]
 
  fit <- tryCatch(
    cfa(cfa_model,
        data = d,
        ordered = items_all,
        estimator = "WLSMV",
        std.lv = TRUE),
    error = function(e) NULL
  )
 
  if (is.null(fit)) return(NA)
 
  if (!lavInspect(fit, "converged")) return(NA)
 
  return(omega_from_fit(fit))
}

# --------------------------
# POINT ESTIMATE
# --------------------------
fit_point <- cfa(cfa_model,
                 data = datos,
                 ordered = items_all,
                 estimator = "WLSMV",
                 std.lv = TRUE)

omega_point <- omega_from_fit(fit_point)

cat("\nOmega:", round(omega_point, 3))

# --------------------------
# BOOTSTRAP
# --------------------------
set.seed(123)

boot_res <- boot(
  data = datos,
  statistic = boot_fun,
  R = 500,
  parallel = "multicore"
)

# clean NAs
boot_vals <- boot_res$t[!is.na(boot_res$t)]

# CI
ci <- quantile(boot_vals, c(0.025, 0.975))

cat("\n\nOmega with 95% CI:\n")
cat(sprintf("%.3f [%.3f, %.3f]", omega_point, ci[1], ci[2]))


SPSS code (version 29.0.0):

* Note: all paths used are absolute; in case of interest in repeat the analyses, the data path must be changed

### KRUSKAL WALLIS ANALYSIS CODE

NPTESTS
  /INDEPENDENT TEST (INTEGRACIÓN SATISFACCIÓN_VITAL BIENESTAR PHQ_TOTAL) GROUP (GASAsinmissing & COREsinmissing & ICDsinmissing)
  /MISSING SCOPE=ANALYSIS USERMISSING=EXCLUDE
  /CRITERIA ALPHA=0.05  CILEVEL=95.

CROSSTABS
  /TABLES=GASAsinmissing & COREsinmissing & ICDsinmissing BY integración10percent satisfaction10percent bienestar10percent depresion10percent
  /FORMAT=AVALUE TABLES
  /STATISTICS=CHISQ PHI
  /CELLS=COUNT EXPECTED ASRESID
  /COUNT ROUND CELL.